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ABSTRACT 



In this article we present results from three on-going 
projects related to the formation of protoplanets in 
protostellar discs. We present the results of simula- 
tions that model the interaction between embedded 
protoplanets and disc models undergoing MHD tur- 
bulence. We review the similarities and differences 
that arise when the disc is turbulent as opposed to 
laminar (but viscous), and present the first results 
of simulations that examine the tidal interaction be- 
tween low mass protoplanets and turbulent discs. 
We describe the results of simulations of Jovian mass 
protoplanets forming in circumbinary discs, and dis- 
cuss the range of possible outcomes that arise in hy- 
drodynamic simulations. 

Finally, we report on some preliminary simulations 
of three protoplanets of Jovian mass that form ap- 
proximately coevally within a protostellar disc. We 
describe the conditions under which such a system 
can form a stable three planet resonance. 
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1. INTRODUCTION 



The on-going discovery of extrasolar planets has 
reinvigorated efforts to understand the formation 
and evolution of planetary systems (Mayor & Queloz 
1995; Marcy & Butler 1996; Marcy, Cochran, & 
Mayor 1999; Vogt et al. 2002; Santos et al. 2003). As 
observations are carried out over longer time scales, 
and with increasing sensitivity, the physical prop- 
erties of the observed planetary systems are set to 
diversify significantly. At the present time, however, 
all of the known planets are Jovian-like gas giants. 
For this reason, much of the theoretical research cur- 
rently underway is examining the various stages of 
formation and evolution of giant protoplanets. 

The most widely accepted theory of how giant plan- 
ets form, the so-called core instability model, sug- 
gests that a multi stage process operates. The solid 



component of protostellar discs gradually coagulates 
to form a solid core of around 15 Earth masses, onto 
which a gaseous envelope accretes (e.g. Bodenheimer 
& Pollack 1986; Pollack et al. 1996). An alterna- 
tive model suggests that giant protoplanets form via 
gravitational instability of a protostellar disc (e.g. 
Boss 2001). In either scenario, interaction between 
the forming protoplanet and the protostellar disc is 
likely to significantly affect the evolution, leading to 
orbital migration. 

In the standard picture of disc-planet interactions, 
a protoplanet exerts torques on a protostellar disc 
through the excitation of spiral density waves at 
Lindblad resonances (e.g. Goldreich & Tremaine 
1979). These waves carry an associated angular mo- 
mentum flux which is deposited in the disc where the 
waves are damped. This process results in a nega- 
tive torque acting on the protoplanet from the outer 
disc and a positive torque acting on it from the disc 
interior to its orbit. For most disc models, the outer 
disc torque is dominant, leading to inward migration 
(Ward 1997). 

A sufficiently massive protoplanet can open up an 
annular gap in a viscous disc centred on its orbital ra- 
dius (Papaloizou & Lin 1984). For typical protostel- 
lar disc models the protoplanet needs to be approx- 
imately a Jovian mass for gap formation to occur. 
Recent simulations (Bryden et al. 1999; Kley 1999; 
Lubow, Seibert, & Artymowicz 1999; D'Angelo, Hen- 
ning, & Kley 2002) examined the formation of gaps 
by giant protoplanets, and also estimated the max- 
imal gas accretion rate onto them. The orbital evo- 
lution of a Jovian mass protoplanet embedded in a 
standard laminar viscous protostellar disc model was 
studied by Nelson et al. (2000). They found that 
that gap formation and accretion of the inner disc 
by the central mass led to the formation of a low 
density inner cavity in which the planet orbits. In- 
teraction with the outer disc resulted in inward type 
II migration on a time scale of a few x 10^ yr. Gas 
accretion during migration allows the protoplanet to 
grow to ~ 3-4 Jupiter masses. The disc models in 
these studies all adopted an anomalous disc viscosity 
modelled through the Navier-Stokes equations with- 
out consideration of its origin. 



The most likely origin of the viscosity is through 
MHD turbulence resulting from the magnetorota- 
tional instability (MRI) (Balbus & Hawley 1991) 
and it has recently become possible through improve- 
ments in computational resources to simulate discs in 
which this underlying mechanism responsible for an- 
gular momentum transport is explicitly calculated. 
This is necessary because the turbulent fluctuations 
do not necessarily result in transport phenomena 
that can be modelled with the Navier-Stokes equa- 
tion. 

To this end Papaloizou & Nelson (2003) and Nelson 
& Papaloizou (2003a) developed models of turbulent 
protostellar accretion discs and considered the inter- 
action with a giant protoplanet of 5 Jupiter masses. 
The large mass was chosen to increase the scale of the 
interaction so reducing the computational resources 
required. This protoplanet was massive enough to 
maintain a deep gap separating the inner and outer 
disc and exert torques characteristic of type II mi- 
gration. More recent work has expanded the range 
of protoplanet masses examined (Papaloizou, Nelson, 
& Snellgrove 2003; Nelson & Papaloizou 2003b). We 
discuss some of the main points of this work in later 
sections. 

The majority of the planets so far detected orbit 
around single solar-type stars, but there have been 
also been detections in binary systems [e.g. 7 Cephei 
(Hatzes et al. 2003), 16 Cygni B (Cochran et al. 
1997)]. Most field stars appear to be members of bi- 
nary systems (Duquennoy & Mayor 1991). For the 
longer period systems planets may orbit around one 
member of the binary. In the shorter period systems 
they could orbit stably around both stars (i.e. cir- 
cumbinary planets). 

The majority of T Tauri stars, whose discs are 
thought to be the sites of planet formation, also 
appear to be in binary or multiple systems, (Ghez, 
Neugebauer & Matthews 1993; Leinert et al. 1993; 
Mathieu et al. 2000). Most have sufficiently large 
separations that it is expected that each component 
will have its own circumstellar disc. For shorter pe- 
riod systems, however, one expects the existence of 
a circumbinary disc, a number of which have been 
observed (e.g. DQ Tau, AK Sco, UZ Tau, GW Ori, 
GG Tau). 

The confirmed existence of planets in binary sys- 
tems, combined with the fact that binary systems 
appear to be common, and to be present during the 
T Tauri phase, means that it is of interest to ex- 
plore how stellar multiplicity affects planet forma- 
tion, and post-formation planetary orbital evolution, 
including formation in circumbinary discs. Previous 
work examined the stability of planetary orbits in 
binary systems using N body simulations (Dvorak 
1986; Holman & Wiegert 1999). This work showed 
that there is a critical ratio of planetary to binary 
semimajor axis for stability, depending on the binary 
mass ratio, qun, and eccentricity eun- A recent pa- 
per (Quintana et al. 2002) explored the late stages of 
terrestrial planet formation in the a Centauri system. 
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Figure 1. This figure shows midplanet density distri- 
bution of 5 Jupiter mass protoplanet in a turbulent 
disc 



This work concluded that the binary companion can 
help speed up planetary accumumlation by stirring 
up the planetary embryos, thus increasing the col- 
lision rate. Recent work by Kley & Burkert (2000) 
examined the effect that an external binary compan- 
ion can have on the migration and mass accretion of 
a giant protoplanet forming in a circumstellar disc. 
They found that for sufficiently close companions, 
both the mass accretion rate and the orbital migra- 
tion rate could be increased above that expected for 
protoplanets forming around single stars. 

The question of how protoplanet evolution is affected 
in circumbinary discs has been examined recently by 
Nelson (2003a) . This work explored the evolution of 
Jovian mass protoplanets forming in circumbinary 
discs using hydrodynamic simulations of a binary 
star plus protoplanet system interacting with a vis- 
cous protostellar disc. The models apply primarily to 
binaries with orbital periods of ~ 1 yr, and semima- 
jor axes of ~ 1 AU that have protoplanets forming 
at a radius of a few AU in the circumbinary disc, al- 
though the results can be scaled to apply to different 
parameters. It is well known that a giant protoplanet 
embedded in a disc around a single star undergoes in- 
ward migration driven by the viscous evolution of the 
disc (e.g. Nelson et al. 2000). The work by Nelson 
(2003a) examines how this process is affected when 
the central star is replaced by a close binary system, 
and delineates the various modes of behaviour that 
arise depending on the properties of the system (e.g. 
binary mass ratio qhim binary eccentricity ei,im etc). 
We present some of the main results of this work in 
later sections. 

A number of the known extrasolar planets exist in 

multi-planet systems. Three of these systems con- 
tain a pair of planets that are in mean motion reso- 
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Figure 2. This figure shows a close-up of the mid- 
plane density distribution for a 5 Jupiter mass pro- 
toplanet in a turbulent disc 

nance (GJ876, HD82943, 55 Cancri). The most likely 
explanation for these resonant systems is approxi- 
mately coeval formation of the planet pair, followed 
by disc-induced differential migration leading to res- 
onant capture (e.g, Snellgrove, Papaloizou, & Nel- 
son 2001; Lee & Peale 2002). In later sections we 
present some preliminary results that examine the 
plausibility of three-planet resonances being discov- 
ered in which the three planets are of approximately 
Jovian mass. The Jovian satellite system displays 
such a configuration with lo, Europa, and Ganymede 
all participating in the Laplace resonance (e.g Peale, 
Cassen, & Reynolds 1979). Here lo and Europa are 
in 2:1 resonance, and Europa and Ganymede are si- 
multaneously in 2:1 resonance, leading to a 4:2:1 re- 
lationship between the mean motion of lo, Europa, 
and Ganymede. Our preliminary calculations in- 
dicate that three planet resonances may indeed be 
established, but that the 4:2:1 relation is unstable. 
However, a situation in which one of the planet pair 
is in a 3:1 resonance, leading to a 6:2:1 or 6:3:1 rela- 
tion between the mean motions, appears to be sta- 
ble, suggesting that such configuration may be found 
among the population of extrasolar planets. 

This article is organised as follows. In section |5] 
we present the results of simulation of high and low 
mass protoplanets interacting with turbulent accre- 
tion discs. In section 01 we describe the results of 
simulations that examine the evolution of giant pro- 
toplanets forming in circumbinary accretion discs. In 
section0]we present some preliminary results of three 
planet systems leading to the formation of three- 
planet resonances. Finally we summarise the results 
presented in this article in section [S] 



2. TURBULENT DISC - PROTOPLANET 
INTERACTIONS 

Most of the previous work examining the interac- 
tion between protostellar disc models and embedded 
protoplanets have used the Navier-Stokes equations 
to simulate laminar disc models with an anomalous 
viscosity coefficient (e.g. Bryden et al. 1999; Kley 



Figure 3. This figure shows magnetic field vectors 
for equivalent region shown in fig .2. 

1999; Nelson et al. 2000). The most likely source 
of disc viscosity, however, is MHD turbulence that 
arises because of the MRI (Balbus & Hawley 1991). 
We present simulations of magnetic, turbulent discs 
interacting with protoplanets of different mass in the 
following sections. 



2.1. Giant Protoplanets 

A detailed discussion of the results presented in 
this section may be found in Nelson & Papaloizou 
(2003a) . The underlying disc model used is described 
in detail in Papaloizou & Nelson (2003). The disc 
was a cylindrical disc model with a locally isothermal 
equation of state such that the effective aspect ratio 
H/R — 0.1 throughout. The inital disc model was 
initiated with a zero-net flux vertical magnetic field 
that varied sinusoidally with radius, and was allowed 
to run until a statistically steady state turbulent disc 
model was obtained. The initial value of the volume 
averaged plasma (3 was < f3 >~ 1000, where (3 is 
the ratio of the thermal gas pressure to the magnetic 
pressure. The final turbulent state generated a vol- 
ume averaged value of the Shakura-Sunyaev stress 
parameter < a >~ 5 x 10^"^. A protoplanet of mass 
Mp ~ 5 Jupiter masses was placed in the disc at a 
radius rp = 2.2. Prior to placing the protoplanet in 
the disc, a partial gap was made in the disc in order 
to avoid the generation of transient features in the 
flow associated with the gap opening process. 

A snapshot of the final state of the disc with the pro- 
toplanet is shown in figure after ^ 100 planetary 
orbits. The protoplanet showed a clear tendency to- 
wards gap formation by deepening and widening the 
initially imposed partial gap during the simulation. 
Broadly speaking, the visual appearance of the disc 
in figure differs significantly from that obtained 
in an equivalent simulation with a laminar, viscous 
disc (e.g. Nelson et al. 2000; Nelson & Papaloizou 
2003a), with the turbulent disc showing a less regu- 
lar and more time dependent structure, in which the 
spiral waves have a more diffused appearence. 

In the region of the wakes the disc and its turbulence 
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Figure 4- This figure shows the azimuthally averaged 
radial density distribution for MHD turbulent disc 
(dotted line), and laminar disc with (solid line) and 
without (dashed line) alpha viscosity 

are strongly perturbed. Magnetic field vectors in the 
disc midplane in the neighbourhood of the planet 
are illustrated in figure |21 along with a correspond- 
ing density plot for the equivalent region in figure El 
An inspection of the magnetic field vectors indicates 
that these tend to line up along the location of the 
wakes but in a somewhat broadened region slightly 
behind the shocks. In this way an ordered structure 
appears to be imposed on the flow and magnetic field 
by the protoplanet. The magnetic stress is largely 
communicated in these ordered regions, leading to 
a significant change in the magnetic contribution to 
the stress parameter a there (Nelson & Papaloizou 
2003a). 

In addition to running the full MHD turbulent disc 
model, we have also run some 2D laminar a disc 
models for the purposes of comparison. Two 2D 
model were run with a = and a = 5 x 10^"^, respec- 
tively. For initial conditions in these 2D models, we 
took the midplane density distribution of the 3D tur- 
bulent model just prior to the addition of the planet, 
switched off the magnetic field, and introduced the 
required value of a in the Navier-Stokes viscosity. 
The implementation of the Navier-Stokes viscosity 
is described in Nelson et al. (2000). 

The azimuthally averaged surface density distribu- 
tion for these models, plus the midplane density dis- 
tribution for the turbulent model, is plotted in fig- 
ure^ Each of the models have been run for a total 
time of ~ 2050 time units, corresponding to ^ 100 
planetary orbits. It is clear that the a — 0.0 2D 
model (dashed line) has the deepest and widest gap, 
as expected. It is also apparent that the turbulent 
model (dotted line) has a deeper and wider gap than 
the 2D a = 5 X lO'^ model, even though a volume 
averaged estimate for the underlying turbulent disc 
model yields an effective a ~ 5 x lO^'^ (see figure 18 
in Papaloizou & Nelson 2003). Thus, the turbulent 
model behaves as if it has a somewhat smaller a than 
reasonable estimates suggest it has. This may arise 
for the following reasons. First, a Navier-Stokes vis- 
cosity with anomalous viscosity coefficient provides a 
source of constantly acting friction in the disc, such 



that it can induce a steady mass flow into the gap 
region. The turbulence, however, does not operate 
as a constant source of friction that generates steady 
inflow velocities. Instead it generates large velocity 
fluctuations that may be much larger than the un- 
derlying inflow velocity arising from the associated 
angular momentum transport. Results presented in 
Papaloizou & Nelson (2003) indicate that a process 
of time averaging the turbulent velocity field is re- 
quired over long time periods before these fluctua- 
tions can be averaged out to reveal the underlying 
mass flow. The disc material in the vicinity of the 
planet experiences periodic high amplitude pertur- 
bations induced by the planet on a time scale much 
shorter than the required averaging time scale, so 
that the disc response is expected to differ from that 
in the case of a disc with Navier-Stokes viscosity. 
A second plausible reason for the apparently lower 
a is that the existence of the magnetic field in the 
turbulent disc allows for field lines to connect across 
the gap region, and to enable angular momentum 
transport across the gap. In this way the magnetic 
field actually helps the planet to maintain the gap 
(Nelson & Papaloizou 2003a). The differences in gap 
structure found here suggest that the accretion rate 
onto a protoplanet in an MHD turbulent disc is likely 
to be less than previous estimates based on laminar 
a disc models indicate [e.g. Bryden et al. (1999); 
Kley (1999); Lubow, Seibert, & Artymowicz (1999); 
Nelson et al. (2000); D'Angelo, Henning, & Kley 
(2002)], but not by a large factor. 

The velocity field within the Hill sphere of the planet 
is shown in the left hand panel of figure|31for the tur- 
bulent disc. An equivalent plot for a laminar disc is 
shown in the right hand panel of this figure. The ma- 
terial that enters the Hill sphere of the protoplanet 
and becomes gravitationally bound to it is normally 
expected to circulate around it by virtue of its angu- 
lar momentum. In the calculation presented here the 
protoplanet does not accrete gas, and is modelled as 
a softened point mass. Consequently it is expected 
that material that enters the Hill sphere will form 
a hydrostatic 'atmosphere' around the planet that 
is supported through pressure and rotation. In fig- 
ure we would expect to see circulation occurring 
in a clockwise fashion for both the MHD turbulent 
disc and the laminar disc. However, it is apparent 
that the circulating pattern has been disrupted in 
the magnetised disc, indicating that magnetic break- 
ing may have been responsible for this (Nelson & 
Papaloizou 2003a). The expected circulation is ob- 
served in the non magnetic run performed at the 
same numerical resolution as shown in the right hand 
panel of figure Magnetic braking of material that 
forms a circumplantary disc inside the Hill sphere 
will have significant implications for the mass accre- 
tion onto the protoplanet. 



2.2. Low mass protoplanets 

In this section we present the results of simulations 
that examine the interaction of turbulent discs with 




Figure 5. This figure shows velocity vectors for material in the planet Hill sphere. The left hand panel shows 
the MHD run, and the right hand panel a laminar disc run. 




Figure 6. This figure shows a snapshot of a turbu- 
lent disc with an embedded q — 3 x 10^^ protoplanet 
located at (x,y)=(-3,0). 



low mass, embedded protoplanets. A detailed discus- 
sion of the results presented in this section is given in 
Papaloizou, Nelson & Snellgrove (2003) and Nelson 
& Papaloizou (2G03b). The underlying disc model 
used in the simulations with low mass protoplanets 
differed from that described in section ITTl The disc 
was a cylindrical disc model with a locally isother- 
mal equation of state. The effective aspect ratio 
took a constant value H/r = 0.07. The disc model 
was initiated with a zero-net flux toroidal magnetic 
field that varied sinusoidally with radius. The initial 
value of the volume averaged plasma [3 parameter 
< fj >~ 30. The disc model was evolved until a sta- 
tistically steady state turbulent disc was obtained. 
The final turbulent state had an associated volume 
averaged stress parameter < a >~ 7 x f0~'^ (Pa- 
paloizou, Nelson, & Snellgrove 2003). 

Once the final turbulent state had been obtained, low 
mass protoplanets were inserted into the disc model. 



These had mass ratios q — mp/M^. — 10 ^, 3 x fO ^, 
and 10"'', respectively. The gravitational potential 
of the protoplanets was softened using a softening 
parameter h — 0.3Hp, where Hp is the disc semi- 
thickness at the position of the protoplanet. The 
primary aim of these simulations was to examine the 
tidal torques exerted on non gap forming protoplan- 
ets by turbulent accretion discs to estimate the cor- 
responding migration rates. 

A snapshot of a simulation with q = 3 x 10^^ is 
shown in figureEl This corresponds to a protoplanet 
of TUp ~ 10 Earth masses orbiting a solar mass star. 
The protoplanet is located at {x,y) = (—3,0), and 
is only just visible in this figure since the density 
fluctuations generated by the turbulence are in fact 
of higher amplitude than the spiral wakes generated 
by the protoplanet. The torque per unit mass ex- 
erted on the protoplanet was calculated as a func- 
tion of time for all simulations, and is shown in fig- 
ure [3 for the case with q = 3 x 10~^. This figure 
clearly shows that the force exerted on the proto- 
planet by the turbulent disc is dominated by high 
amplitude fluctuations. The total torque exerted on 
the protoplanet is seen to oscillate between positive 
and negative values, suggesting that the migration 
in this case is likely to occur as a 'random walk', 
rather than as a monotonic inward drift normally 
associated with type I migration (Ward 1997). Fig- 
ure IHl shows the time evolution of the running time 
average of the torque per unit mass. The straight 
line shows the torque obtained from an equivalent 
simulation with a laminar disc model. The upper 
line shows the running time average of the torque 
due to the inner disc, the lower line the torque due 
to the outer disc, and the middle line the running 
average of the total torque. The running mean of 
the total torque is apparently not converging to a 
well defined value on the time scale of the simula- 
tion, and for a large part of the simulation indicates 
that the protoplanet would migrate outwards on av- 
erage. Treating the turbulent fluctuations as having 
a Gaussian distribution, we can estimate the time 
for the running mean to converge as a function of 
the amplitude of the fluctuations. Such an estimate 
yields a time scale of ~ 100 planetary orbits for con- 
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Figure 7. This figure shows the time evolution of the 
torque per unit mass exerted on q ^ 3 x 10^^ pro- 
toplanet. The torque is clearly dominated by strong 
fluctuation due to the turbulence. 

vergence, longer than we are currently able to evolve 
the simulation (Nelson & Papaloizou 2003b). 

The simulations with q — 10^^ and q — 10~^ showed 
qualitatively similar results to those described above. 
Although we are unable to run simulations for suf- 
ficient length to definitively calculate the migration 
times of protoplanets in turbulent discs, it is clear 
that the picture of migration that emerges differs sig- 
nificantly from that obtained in laminar disc models. 
These results offer the possibility that the rapid mi- 
gration rates of giant protoplanet cores may in fact 
be slowed down or stopped by interaction with the 
density fluctuations in a turbulent disc. 



3. GIANT PLANETS FORMING IN 
CIRCUMBINARY DISCS 

In this section we present the results of simulations 
that examine the evolution of giant protoplanets that 
form in circumbinary discs. A fuller discussion of this 
work is presented in Nelson (2003a). We consider 
the interaction between a coplanar binary and pro- 
toplanet system and a two-dimensional, gaseous, vis- 
cous, circumbinary disc within which it is supposed 
the giant protoplanet formed. We do not consider 
the early evolution of the protoplanet in this work 
or address the formation process itself, but make the 
assumption that a Jovian mass protoplanet is able 
to form and examine the dynamical consequences of 
this. 

Each of the stellar components and the protoplanet 
experience the gravitational force of the other two, as 
well as that due to the disc. The disc is evolved using 
the hydrodynamics code NIRVANA (Ziegler & Yorke 
1997). The planet and binary orbits are evolved us- 
ing a fifth-order Runge-Kutta scheme (Press et al. 



Time Averaged Torques Per Unit Mass (Mp— 10) 




5 10 15 

Time (Orbits) 



Figure 8. This figure shows the time evolution of the 
running time average of the torque per unit mass. 
The upper line corresponds to the inner disc torque, 
the lower line gives the outer disc torque, and the 
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Figure 9. This figure shows surface density contours 
for run in which the planet is ejected by the binary 



1992). The force of the planet on the disc, and of 
the disc on the planet, is softened using a gravita- 
tional softening parameter b = 0.5ap{H/r), where Op 
is the semimajor axis of the planet, and H/r is the 
disc aspect ratio. We assume that the mass of the 
protoplanet is fixed, and so do not allow accretion of 
matter from the disc onto the protoplanet. 
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We adopt a disc model in which the effective aspect 
ratio H/r ~ 0.05, and the Shakura-Sunyaev viscos- 
ity parameter a = 5 x 10^'^ (Shakura & Sunyaev 
1973). The surface density S is initialised to have 
an inner cavity within which the planet and binary 
orbit (Nelson 2003a). Simulations initiated with no 
inner cavity show that one is formed by the action of 
the binary system and planet clearing gaps in their 
local neighbourhood. As the planet migrates in to- 
wards the central binary these gaps join to form a 
single cavity. The disc mass is normalised through 
the choice of Sq such that a standard disc model with 
S(r) — Tior^^^^ throughout would contain about 4 
Jupiter masses interior to the initial planet radius 
(assumed in physical units to be 5 AU). Thus the 
disc mass interior to the initial planet radius would 
be about twice that of a minimum mass solar nebula 
model. 

The total mass of the binary plus protoplanet sys- 
tem is assumed to be 1 Mq. Dimcnsionless units 
are used such that the total mass of the binary sys- 
tem plus planet Mtot = 1 and the gravitational con- 
stant G = 1. The initial binary semimajor axis is 
cibin = 0.4 in all simulations, and the initial planet 
semimajor axis ap = 1.4. The simulations are initi- 
ated with the binary system having an initial eccen- 
tricity, Cbin, and the protoplanet is initially in cir- 
cular orbit. The binary mass ratio qbm = 0.1 for 
all simulations presented in this section, but larger 
values were considered in Nelson (2003a). The unit 
of time quoted in the discussion of the simulation 
results below is the orbital period at i? = 1. 
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Figure 10. This figure shows semimajor axes and 
eccentricities for run in which planet is scattered by 
the binary. 



principle the protoplanet is able to enter higher order 
resonances than 4:1, such as 5:1 or 6:1, since its initial 
location lies beyond these resonance locations. How- 
ever, none of the simulations presented here resulted 
in such a capture. Test calculations indicate that 
capture into higher order resonances requires slower 
planetary migration rates than those that arise in 
these simulations. For significantly faster migration 
rates the planet may pass through the 4:1 resonance 
(Nelson 2003a). 



3.1. Numerical Results 

The results of the simulations can be divided into 
three categories (Mode 1, Mode 2, and Mode 3), 
which are described below, and are most strongly 
correlated with changes in the binary mass ratio, 
qbim and binary eccentricity Cbin- Changes to the 
disc mass and/or protoplanet mass appear to be less 
important. In some runs the planet enters the 4:1 
mean motion resonance with the binary. The associ- 
ated resonant angles in the coplanar case are defined 
by: 



3.1.1. Mode 1 - Planetary Scattering 

A number of simulations resulted in a close encounter 
between the protoplanet and binary system, lead- 
ing to gravitational scattering of the protoplanet to 
larger radii, or into an unbound state. We label this 
mode of evolution as 'Mode 1'. Typically the ini- 
tial scattering causes the eccentricity of the planet 
to grow to values Cp ~ 0.9, and the semimajor axis 
to increase to Op ~ 6 - 8. In runs that were contin- 
ued for significant times after this initial scattering, 
ejection of the planet could occur after subsequent 
close encounters. 



i'l — 4 As — Ap — 3uJs ip2 = 4As — Ap — 3a;p (1) We will illustrate this mode of evolution using a sim- 

= — Ap — 2ujs — ujp ip4 — 4As — Ap — 2a;p — ujs ulation with nip = 3 Jupiter masses and qbin = 0.1. 

A series of snapshots of the simulation are shown in 

where A^, Ap are the mean longitudes of the sec- figure El Mode 1 evolution proceeds as follows. The 

ondary star and protoplanet, respectively, and Wg, Wp protoplanet migrates in towards the central binary 

are the longitudes of pericentre of the secondary and due to interaction with the circumbinary disc, and 

protoplanet, respectively. When in resonance tJjt, or temporarily enters the 4:1 mean motion resonance 

■04 should librate, or all the angles should librate. In with the binary. The migration and eccentricity evo- 





Figure 11. This figure shows the evolution of the 
resonant anglesd defined in equation\^ 

lution is shown in figure IIUI and the resonance an- 
gles are shown in figure 1111 The resonant angle V'a 
librates with low amplitude, indicating that the pro- 
toplanet is strongly locked into the resonance. The 
resonance drives the eccentricity of the protoplanet 
upward, until the protoplanet has a close encounter 
with the secondary star during or close to periapse, 
and is scattered out of the resonance into a high ec- 
centricity orbit with significantly larger semimajor 
axis. We note that being in resonance normally helps 
maintain the stability of two objects orbiting about 
a central mass. However, when one of the objects 
is a star, the large perturbations experienced by the 
planet can cause the resonance to break when the 
eccentricities are significant. Once out of resonance, 
the chances of a close encounter and subsequent scat- 
tering are greatly increased. This provides a method 
of forming 'free-floating planets'. 



3.1.2. Mode 2 - Near-resonant Protoplanet 

A mode of evolution was found in some of the simula- 
tions leading to the protoplanet orbiting stably just 
outside of the 4:1 resonance. We label this mode of 
evolution as 'Mode 2'. Mode 2 evolution is illustrated 
by a simulation for which rup — 1, qhin = 0.1, and 
Sbin = 0.1. The evolution of the orbital elements 
are shown in figure 1121 Here, the protoplanet mi- 
grates inwards and becomes weakly locked into the 
4:1 resonance, with the resonant angle ips librating 
with large amplitude. The resonance becomes un- 
defined and breaks when ep — momentarily dur- 
ing the high amplitude oscillations of ep that accom- 
pany the high amplitude librations of ^3. The pro- 
toplanet undergoes a period of outward migration 
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Figure 12. This figure shows semimajor axes and 
eccentricities for the Mode 2 run described in the text. 

through interaction with the disc by virtue of the 
eccentricity having reattained values of Cp ~ 0.17 
once the resonance is broken. Calculations by Nelson 
(2003b) show that gap-forming protoplanets orbiting 
in tidally truncated discs undergo outward migration 
if they are given eccentricities of this magnitude im- 
pulsively, due to the sign of the torque exerted by the 
disc reversing for large eccentricities. The outward 
migration moves the planet to a safer distance away 
from the binary, helping to avoid instability. 

Once the protoplanet has migrated to just be- 
yond the 4:1 resonance the outward migration halts, 
since its eccentricity reduces slightly, and the planet 
remains there for the duration of the simulation. 
The system achieves a balance between eccentricity 
damping by the disc and eccentricity excitation by 
the binary, maintaining a mean value of Cp ~ 0.12 
(Nelson 2003a). The torque exerted by the disc on 
the protoplanet is significantly weakened by virtue 
of the finite eccentricity (Nelson 2003b), preventing 
the planet from migrating back towards the binary. 
Continuation of this run in the absence of the disc in- 
dicates that the planet remains stable for over 6 x 10^ 
orbits. This is in good agreement with the stability 
criteria obtained by Holman & Wiegert (1999) since 
the protoplanet lies just outside of the zone of insta- 
bility found by their study. 



3.1.3. Mode 3 - Eccentric Discs 

A mode of evolution was found in which the plan- 
etary migration was halted before the protoplanet 
could approach the central binary and reach the 4:1 
resonance. This only occurred when the central bi- 
nary had an initial eccentricity of Chin > 0.2. The mi- 
gration stalls because the circumbinary disc becomes 
eccentric. We label this mode of evolution as 'Mode 
3'. We illustrate this mode of evolution by present- 
ing the results of a simulation with rUp = 1 Jupiter 
mass, qun = 0.1, and Cbin = 0.2. Figure [T^ shows 
snapshots of the surface density at different times 
during the simulation, with the disc becoming notice- 
ably eccentric. Interaction between the protoplanet 
and the eccentric disc leads to a dramatic reduction 
or even reversal of the time-averaged torque driv- 
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Figure 13. This figure shows contours of surface den- 
sity for the Mode 3 run described in the text. 



Figure 14- This figure shows semimajor axis and 
eccentricity for Mode 3 run described in text. 

ing the migration. This is because the disc-planet 
interaction becomes dominated by the m = 1 sur- 
face density perturbation in the disc rather than by 
the usual interaction at Lindblad resonances in the 
disc. Figure IT^ shows the evolution of the semima- 
jor axis and eccentricity of the planet, and shows 
that the migration stalls. Simulations of this type 
can be run for many thousands of planetary orbits 
without any significant net inward migration occur- 
ring. Such systems are likely to be stable long after 
the circumbinary disc has dispersed, since the plan- 
ets remain in the region of stability defined by the 
work of Holman & Wiegert (1999), and are probably 
the best candidates for finding stable circumbinary 
extrasolar planets. Interestingly, spectroscopic bi- 
nary systems with significant eccentricity are signif- 
icantly more numerous than those with lower eccen- 
tricities (e.g. Duquennoy & Mayor 1991; Mathieu et 
al. 2000), suggesting that circumbinary planets may 
be common if planets are able to form in circumbi- 
nary discs. 



4. THREE PLANET RESONANT SYSTEMS 

The are a number of multiplanet systems among the 
known extrasolar planets. In at least two, and possi- 
bly three, of these multiplanet systems, two planets 
appear to be in mean motion resonance. The plan- 
ets in the systems GJ876 and HD82943 appear to be 
in 2:1 resonances, while two of the planets in the 55 
Cancri system appear to be in 3:1 resonance. The 
existence of these mean motion resonances can be 
explained by a model in which two planets form in 
the disc approximately coevally, and are driven into 
resonance by differential inward migration due to in- 
teraction with the protostellar disc (e.g. Snellgrove, 
Papaloizou, & Nelson 2001; Lee & Peale 2002; Nel- 
son & Papaloizou 2002). 

In this section we present some preliminary calcu- 
lations that address the question of whether three- 
planet resonances may be formed, and remain sta- 
ble, by a process which is analogous to the above 
scenario when each of the protoplanets is of Jovian 
mass. Namely, approximately coeval formation of 
three planets which are then driven into a three- 
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Figure 15. This figure shows contour plots for the 
run designed to test scenario (i) described in text. 

planet resonance through interaction with the proto- 
planetary disc. We examine this question by means 
of 2~D numerical simulations very similar to those 
described in section fusing the same code and basic 
disc models. 

We have considered a number of different scenarios, 
which we describe below: 

(i) Coeval formation of three planets, with the outer 
two forming a 2:1 resonance, followed by inward mi- 
gration of the resonant pair until the middle and in- 
ner planet themselves become locked in a 2:1 reso- 
nance. We refer to this as a 4:2:1 resonance. 
(m) Formation of two planets which differentially mi- 
grate and lock into 2:1 resonance, followed by forma- 
tion of a third planet at larger radius that migrates in 
and locks into 2:1 resonance with the middle planet. 
This is another example of a 4:2:1 resonance, 
(m) A scenario very similar to scenario (z), except 
the middle and inner planet form a 3:1 resonance. 
We call this a 6:2:1 resonance. 

{iv) A scenario similar to scenario {ii), except the 
outer and middle planets form a 3:1 resonance. We 
call this a 6:3:1 resonance. 



4.1. 4:2:1 Resonances 

Both scenrios ( i) and [ii) above could potentially lead 
to the formation of a 4:2:1 resonance. In this sub- 
section we present the results of a simulation that 
examines scenario («). This simulation consisted of 
three Jovian mass protoplanets being inserted into 
a protostellar disc model on fixed circular orbits for 
a period of ~ 200 orbits of the middle planet until 
partial gap clearing had occurred. The protoplanets 
were then allowed to evolve under the gravitational 
influence of the disc. The evolution is shown in fig- 
ure 1151 which shows snapshots of the evolution at 



Semimajor a^is versus time 




□ 1000 2000 5000 4000 

time (orbits) 



Figure 16. This figure shows semimajor axes and 
eccentricities for three planets in run designed to test 
scenario (i) described in text. 



four separate times. The orbital evolution is shown 
in figure 1161 The outer planet migrates inward the 
most rapidly and enters the 2:1 resonance with the 
middle planet. This resonant pair then migrate in 
towards the innermost planet. As they approach the 
orbital radius at which the middle planet becomes 
locked into 2:1 resonance with the innermost planet, 
there is a strong gravitational interaction that leads 
to scattering of the planets. Figure El shows that 
the middle and outer planet interchange during this 
encounter leading to the middle planet being flung 
out to larger radii in the disc. All simulations that 
examine scenario (i) result in similar behaviour, with 
no 4:2:1 resonance ever forming. However, interest- 
ing results are obtained which show that the kind 
of interaction described here may help some planets 
survive migration in a disc when they would other- 
wise migrate into the central star by being ejected to 
larger radii. The evolution after ejection is likely to 
involve gap formation and slow inward type II mi- 
gration. 



Simulations that test scenario (ii) suggest that 4:2:1 
resonances can be established, but that they are 
never stable for more than a few hundred orbits This 
is primarily because the middle planet has its eccen- 
tricity pumped up to large values by virtue of under- 
going resonant interaction with both the inner most 
and outer most planets. These will be described in 
a forthcoming paper (Nelson & Papaloizou 2003c). 




Figure 17. This figure shows contour plots for run 

designed to examine scenario (Hi) described in text. Figure 19. This figure shows the resonant angles 

(upper two panels) associated with 2:1 resonance for 
outer two planets in scenario (Hi) run described in 
text. The lowe.tt nnn.el shows the difference in the 
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Figure 18. This figure shows semimajor axes and 
eccentricities for three planets in run designed to ex- 
amine scenario (Hi) described in text. 



Figure 20. This figure shows the resonant angles 
(upper three panels) associated with 3:1 resonance 
for middle and inner planet in scenario (Hi) run de- 
scribed in text. The lowest panel shows the difference 
in the longitude of pericentre for the middle and in- 
ner planet in this run. 
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4.1.1. 6:2:1 and 6:3:1 Resonances 

Simulations have been performed to examine the for- 
mation of 6:2:1 and 6:3:1 resonances. These indicate 
that such resonances may be formed quite easily, and 
remain stable over long time periods. In this section 
we present one simulation that examined the plausi- 
bility of scenario (iii). 

The initial setup was similar to that described in 
section 14.11 except that the inner two planets were 
more widely separated. The outer most planet mi- 
grated inward the most rapidly and locked into 2:1 
resonance with the middle planet. The resonant pair 
then migrated in towards the inner most planet until 
the middle and inner planet became locked into 3:1 
resonance. FigurelTTIshows snapshots of the disc and 
planet evolution. The orbital evolution is shown in 
figure lTHl which shows the converging orbits, and the 
resulting slow-down in migration as the three planets 
become resonantly locked leading to an effective in- 
crease in the inertia of the system that is being driven 
inward by the disc. The resonant angles for the outer 
two planets are shown in figure El showing the es- 
tablishment of the 2:1 resonance, and the resonance 
angles for the inner two planets are shown in figure !^ 
showing the establishment of the 3:1 resonance. Also 
shown in the lowest panels of figures and EUl are 
the differences in the longitudes of pericentre for each 
of the planet pair. Figure EOl shows that the middle 
and inner planet becomes apsidally locked when the 
outer two planets enter the 2:1 resonance, since the 
eccentricity driving by the resonance acts as an ef- 
fective 'tuning' mechanism for the establishment of 
a secular resonance between the inner two planets. 
We note that such a secular resonance is observed in 
the Upsilon Andromeda system. 



5. SUMMARY 

In this article we have presented the results from 
three distinct projects. The first examined the in- 
teraction between embedded protoplanets and tur- 
bulent, magnetised discs. Broadly speaking we find 
rather similar behaviour for high mass, gap form- 
ing planets when compared with previous work on 
laminar but viscous discs. However, low mass proto- 
planets behave rather differently, and appear to un- 
der go 'stochastic migration' due to interaction with 
the background turbulence. The result appears to 
be that low mass protoplanets will undergo a randon 
walk in turbulent discs, instead of monotonic inward 
migration. 

We presented the results of simulations that exam- 
ined the formation of giant protoplanets in circumbi- 
nary discs. We find that under a wide range of con- 
ditions stable circumbinary planets may be main- 
tained. In particular, if the binary system is signifi- 
cantly eccentric, then disc induced inward migration 
of the protoplanet does not occur. 
Finally, we presented some preliminary calculations 
of three protoplanets forming in a disc, and examined 
the conditions under which three-planet resonances 



could be established and maintained. We found that 
4:2:1 configurations could be formed, but quickly be- 
came unstable. However, 6:2:1 and 6:3:1 configura- 
tions could form and remain stable over very long 
periods. 
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